Nature and number of distinct phases in the random field Ising model 
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We investigate the phase structure of the random-field Ising model with a bimodal random field 
distribution. Our aim is to test for the possibility of an equilibrium spin-glass phase, and for replica 
symmetry breaking (RSB) within such a phase. We study a low-temperature region where the spin- 
glass phase is thought to occur, but which has received little numerical study to date. We use the 
exchange Monte-Carlo technique to acquire equilibrium information about the model, in particular 
the P(q) distribution and the spectrum of eigenvalues of the spin-spin correlation matrix (which tests 
for the presence of RSB). Our studies span the range in parameter space from the ferromagnetic to 
the paramagnetic phase. We find however no convincing evidence for any equilibrium glass phase, 
with or without RSB, between these two phases. Instead we find clear evidence (principally from 
the P(q) distribution) that there are only two phases at this low temperature, with a discontinuity 
in the magnetization at the transition like that seen at other temperatures. 



I. INTRODUCTION 

The random-field Ising model (RFIM) is one of the 
simplest problems that one can define theoretically in 
which disorder is present and plays a significant role 
in the physics. This problem gains further interest 
since a related random system which may be realized 
experimentally — namely, the site-diluted antiferromag- 
net in a uniform field — may be mapped onto the RFIM. 
The RFIM has been subjected to considerable study, 
both experimental and theoretical, for over 25 years [Q. 
Much effort has been focussed on the nature of the phase 
diagram. It is known rigorously that, for low temper- 
atures and low random fields, the ferromagnetic (FM) 
phase found in the pure system survives for d > 3 (where 
d is the space dimension) . However most other questions 
that one can ask about this system for d = 3 (hence be- 
yond the range of validity of mean-field theory) remain 
unanswered. For example, is the phase transition first- 
order or continuous? Is there a tricritical point, separat- 
ing a critical line from a line of first-order transitions? 
Does the answer to these questions, and/or the univer- 
sality class of the critical line, depend on the distribution 
of random fields? And finally: is there a third, spin-glass 
(SG) phase, in addition to (and likely separating) the FM 
and paramagnetic (PM) phases? 

In this work we focus on this last question. Early 
numerical studies using mean-field theory defined 
a noncquilibrium "domain state" , characterized by var- 
ious forms of irreversible behavior, of a similar nature 
to the irreversible behavior seen experimentally. This 
irreversibility region, found between the PM and FM 
phases, gives (at most) hints of the presence of an equi- 



librium glassy phase. Stronger arguments for a glassy 
phase were made by de Almeida and Bruinsma jf|, who 
showed that, although mean- field theory in strictly in- 
finite dimension gives replica symmetry in large di- 
mensions there is a stable glassy phase with replica sym- 
metry breaking (RSB). This result is supported by the 
work of Mezard and Young Q| and of Mezard and Monas- 
son who showed that there is a SG phase with RSB 
in the limit of a large number of spin components, and 
by Guagnelli et al who showed the presence of mul- 
tiple solutions to the mean- field equations in d — 3, at 
a temperature well above the FM ordering temperature. 
Subsequently, Brezin and de Dominicis |l0| ] developed a 
renormalization-group approach for a replicated theory, 
and argued for RSB on the basis of a loss of stability 
of the fixed point. Finally, we mention exact numeri- 
cal zero-temperature studies |PIf which show the appear- 
ance of extensive entropy in the ground state at the point 
where FM order disappears. 

Thus there are good reasons to consider the possibil- 
ity of a SG phase for the RFIM, located between the 
FM and PM phases, and likely (see ||) more easily de- 
tected at low T and large random field. It is also of great 
interest to test for the presence of RSB in this phase. 
To date, conclusive evidence for RSB is confined to the- 
oretical models in infinite dimension; to confirm RSB 
in a three-dimensional problem — which can be studied 
experimentally — would be an important step. 

In earlier work Jl^,[l3| we have developed a method 
for detecting RSB in finite-size numerical studies. The 
method relies on examining the eigenvalues of the spin- 
spin correlation matrix. We showed that, in the limit of 
large number of spins N, there are multiple O(N) eigen- 
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values of this matrix, if |T^] and only if JlJ] there is RSB. 
This was clearly confirmed by Monte-Carlo (MC) studies 
of the infinite-ranged Sherrington-Kirkpatrick problem, 
which is known to have RSB. Unfortunately, this method 
gave no evidence for RSB in a finite- (4)-dimensional ver- 
sion of the same problem. 

We apply the same method to the 3-dimensional RFIM 
in this work. We note that previous numerical studies 
fuj of phase transitions of the RFIM are almost entirely 
confined to two regions: either relatively high tempera- 
tures and low random fields, or zero temperature. How- 
ever, based on the location of the irreversibility regions 
in phase space, and on the results of Ref. ||, one is most 
likely to find a glassy phase by looking in the high-field, 
low-T region; this we do here. 

Our results, in short, are as follows. We find some, 
rather ambiguous, evidence for RSB in the behavior of 
the spin-spin eigenvalues, over a small range of field, at 
fixed, low temperature. On either side of this region we 
see clearly the PM and FM phases. Thus this narrow 
region is either the sought-for glassy phase, or simply 
the finite-size effects of the phase transition itself. Ex- 
amination of the order parameter function P(q) enables 
us to rather conclusively rule out the former explanation 
in favor of the latter: a FM/PM phase boundary. Fi- 
nally, we see at the transition the discontinuity (or near- 
discontinuity |l5|,[l6| ) in magnetization that has been seen 
in many other studies, both at zero JTt|— |l9|] and at higher 
p0|-p2[ temperature. 

We organize the rest of this paper as follows. In Sec. 
II we briefly explain the general technique of correlation 
matrix spectral analysis introduced in our previous work 
| Q |. In Sec. Ill we present the procedures and results of 
the numerical analysis of the bi-modal RFIM model in 
a cross-section of the phase diagram. Finally in Sec. IV 
we present our conclusions. 

II. CORRELATION MATRIX SPECTRAL 
ANALYSIS 

In considering frustrated and disordered systems such 
as spin glasses, one is forced to broaden the concept of 
'ordering' in describing the frozen phase. The key point, 
as recognized long ago by Edwards and Anderson [p3| , 
is simply the freezing itself: long-time averages become 
nonzero. One can detect this freezing through correla- 
tion functions: if the freezing is long-ranged (i.e., truly 
a distinct phase), then so are the correlations. Viewing 
the correlation function as a matrix, one finds that order- 
ing (freezing) appears in the spectrum of the correlation 
matrix as an extensive (which we will also call 'large') 
eigenvalue p4| ]. Thus examination of this spectrum al- 
lows one to detect ordering/freezing, while completely 
avoiding any need to guess the nature of the ordering 



(i.e., the eigenfunction corresponding to the large eigen- 
value). 

For Ising systems, a suitable correlation function is 

Cij = (SiSj) ; 

this is a real symmetric positive semi-definite matrix with 
trace equal to the system volume N. In the frozen phase 
at least one of the eigenvalues of Cy is extensive, i.e., of 
O(N). For a simple Ising ferromagnet without quenched 
disorder, freezing leads to the partitioning of spin config- 
uration space into two disjoint regions (+ and — ); each 
is frozen, and each has the same thermodynamic weight 
since they are related by a symmetry of the Hamilto- 
nian, namely spin inversion. Since Cij is invariant un- 
der this symmetry, the two (symmetry breaking) frozen 
states give only a single large eigenvalue. 

In the case of replica symmetry breaking or RSB, there 
are multiple frozen states which are not related by any 
symmetry of the Hamiltonian. We have shown that, 
when there is RSB, there must be more than one ex- 
tensive eigenvalue in the spectrum of Cy — assuming of 
course that the thermal average is taken over the entire 
configuration space. We have also shown jl2) that, in 
the absence of RSB, there can only be one such exten- 
sive eigenvalue. Given these results, one can obtain a 
definitive answer to the question of whether or not a sys- 
tem exhibits RSB, by determining whether there is or is 
not more than one large eigenvalue in the spectrum. 

This latter question is well defined, but is vulnerable to 
finite-size uncertainties (as are other criteria such as the 
overlap distribution P(q)). Here such uncertainties arise 
because the definition of an extensive eigenvalue depends 
on taking the large- N limit. Thus, in numerical studies 
of finite systems, one must look for convincing evidence 
that one has indeed found the asymptotic behavior, at 
least of the two largest eigenvalues of C, as N increases. 
In particular, one must determine if the second largest 
eigenvalue A2 grows linearly with N at large N. In the 
absence of RSB, A 2 can grow with N as A 2 ~ N 1 - 6 

with S > 0. RSB, in contrast, requires that 5 be strictly 
zero. 

Let us now discuss these ideas as applied specifically 
to the RFIM. In the paramagnetic phase (for nonzero 
values of the random field) there is freezing of the spins, 
which (on average) simply follow the random field. Thus 
the frozen configuration is unique, giving a single large 
eigenvalue. In the ferromagnetic phase there is a + and 
a — state. However, due to the random field, the two 
are no longer strictly related by symmetry; nor — due to 
a net field of 0(N 1 ^ 2 ) — do they have equal thermody- 
namic weight. (These two states are termed 'similar but 
incongruent' by Huse and Fisher |25|.) Hence we again 
expect a single extensive eigenvalue in the FM phase. 

In this work we will look for signs of a spin-glass phase 
in a low-temperature region between the PM and FM 
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phases. If the SG phase has RSB — and we are able to 
study a range of N beyond the threshold |l2j for observ- 
ing RSB — then we should find multiple large eigenvalues 
of C. Note that, if X2/N decays over the entire range 
of N studied, and reaches a very small value within that 
range, then that small value offers a rough upper bound 
for the thermodynamic weight which may be assigned to 
any putative RSB state. Hence there may be RSB in 
such a case; but it would be a 'weak' RSB, where only 
one state has large thermodynamic weight. We also do 
not rule out the possibility of a SG phase without RSB, 
i.e., with a single large eigenvalue. In this case one needs 
other measures (such as qualitative changes in the scal- 
ing of A2/-/V, in P(q), or in the spin-glass susceptibility 
Xsg) to distinguish the SG phase from the other two. 

III. NUMERICAL PROCEDURES AND RESULTS 
A. The numerical method 

The RFIM Hamiltonian is given by 

H=-J^2S i S j -^2h i Si (1) 

with (ij) indicating the sum over nearest neighbors, J 
being the strength of the ferromagnetic coupling, and hi 
being a quenched random variable. In this work we will 
take J = 1 , so that temperatures and fields are measured 
in units of J. The two common probability distributions 
for hi are a Gaussian distribution, with variance A, or 
the bi-modal distribution 

P(hi) = ±[6(A - hi) + 6(A + hi)] . 

As we noted in the Introduction, there is some uncer- 
tainty as to which features of the phase diagram are sen- 
sitive to the choice between these two distributions. For 
example, Aharony p6[ | predicts a tricritical point for the 
bimodal distribution, but not for the Gaussian. We will 
work with the bimodal distribution, principally because 
it enables faster Monte-Carlo simulation. Also, we know 
of no reason to expect that a SG phase is more or less 
likely to occur for one or the other distribution. 

As is the case in all glassy systems, the system sizes of 
Monte Carlo simulations are severely restricted by long 
relaxation or equilibration times. The typically long time 
that it takes to jump across large energy barriers in con- 
figuration space can be greatly reduced by the paral- 
lel tempering or exchange Monte Carlo technique [p7| . 
This technique consists of running simultaneously mul- 
tiple replicas of the system at different temperatures, 
while allowing the swapping of spin configurations be- 
tween adjacent temperatures. Such swaps are attempted 



every N s Monte Carlo steps. By choosing the probabil- 
ity of swapping correctly — meaning that the probability 
of a given spin configuration being accepted from an- 
other replica at a different temperature is constrained 
to be consistent with the Boltzmann distribution at the 
original temperature — all copies of the system remain in 
equilibrium under the swapping process. This effectively 
enables the jumping of free-energy barriers, since differ- 
ent replicas are likely to be in different local minima of 
the free energy. 

Another important aspect of Monte Carlo simulations 
in quenched-disorder systems is a set of robust equilibra- 
tion criteria, needed to obtain reliable numerical data. 
The parallel tempering technique does not allow one to 
use the equilibrium criteria used in Reference (^8|, in 
which two methods of calculating the disconnected spin 
glass susceptibility approach the equilibrium value from 
different directions. In the case of spin glasses with Gaus- 
sian random bonds, a different technique was used to 
determine the Monte Carlo time needed for equilibrated 
results (^9|. This criterion cannot be used in the RFIM 
problem, but another one, which can be derived in a sim- 
ilar fashion, can be used for Gaussian distributed random 
fields. In the spirit of reference p9| , it is straightforward 
to show that 

1 

where q E A = (1/N) SJ^tW is tlic Edwards- 
Anderson order parameter, and (• • -}t and [• • -) av indi- 
cate thermal and disorder average respectively. One can 
calculate qEA numerically using the overlap of two un- 
corrected replicas (a and b) at the same temperature, 

^i^EEt'wsfw (3) 
Nt ° i t =i 

or by using expression (|J) to obtain what we will call 

(2) 

q EA . These two quantities approach the thermodynamic 
equilibrium value qEA from opposite directions as t , the 
number of Monte Carlo steps used, increases. Hence one 
can use their convergence to a common value as a crite- 
rion for the equilibration time. 

We have not found an expression analogous to (||) for 
the case of a bi-modal distribution of fields. We note how- 
ever that the time evolution of 9^4 (^o) for the bimodal 
case follows that for the Gaussian distribution (using the 
same parameters) very closely. Therefore, we can take 
the equilibrating value of to obtained for the Gaussian 
distribution, and use it for the bi-modal case as a mini- 
mum Monte Carlo time needed to reach equilibrium for 
each disorder configuration. We retain the bi-modal dis- 
tribution for the bulk of our MC studies as it allows for 
a faster algorithm. 

In Fig. ^ we show the behavior of the two versions of 
qEA, calculated with a Gaussian distribution for T/J — 
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2.2, A/ J = 0.4, and L = 5. We also show these two quan- 

tities for the bi-modal model. It is clear that q E ' A is not 

equivalent to q EA m this case - However, as noted above, 

both q EA and q EA for the bi-modal model plateau at the 

same point where q EA ano - 9ea mcc t f° r the Gaussian 
model. 

As criteria for equilibrium we then have several require- 
ments. First, the equilibration time must be at least 
twice the time needed to establish convergence of q EA 

(2) 

and q EA for the Gaussian case. Second, we require that 

two distinct ways of computing Xsg = F £y ('% i ^.7')t> 
one using two replicas at the same temperature and the 
other a direct Monte Carlo method, agree within statis- 
tical error 

Our third criterion has to do with the exchange Monte 
Carlo method. In order to ensure that equilibrium holds 
for all the various temperatures simulated in parallel, we 
require that the probability histogram for visiting the dif- 
ferent temperatures is fiat, and also that the acceptance 
ratio for swapping among the configurations at different 
temperatures is at least 0.3. In order to avoid the bot- 
tlenecks mentioned in reference (^9| with regard to the 
acceptance ratios among the different temperatures, we 
follow a procedure introduced by Hukushima Here 
we perform a quick Monte Carlo simulation (only a few 
disorder realizations) with an equidistant set of /3's to 
obtain an approximation for the average energy of the 
system. In order to obtain homogeneous acceptance ra- 
tios one performs the following mapping: 

Pn(t + l) = ^\f3 n (t)+g((3 n (t))) , (4) 

where 

-A,+iE03„ +1 ) - E(A,)(/3„_, - . (5) 

Here the values of E((3i) are obtained by extrapolation 
from the ones obtained in the initial Monte Carlo runs. 
Using this procedure we obtain a set of temperatures 
(/3's) that yield homogeneous acceptance probabilities 
when a swap between two adjacent temperatures is at- 
tempted. We can then, if necessary, increase the number 
of /3's and repeat the procedure, until the uniform swap- 
ping probability reaches our minimum goal of at least 
0.3. In Fig. |^ we show the swapping acceptance ratios 
after zero, one, and two iterations of this procedure, for 
L = 11, T m i n /J = 1.5, and A/J — 0.6. It is clear from 
the figure that the bottleneck at T/ J = 4 is removed by 
the iterative procedure. 

Our final criterion is to require that the maximum tem- 
perature used in the parallel tempering process shows 
a 'melted' P{q); that is, P(q) at the smallest j3 must 
have an approximate Gaussian shape centered at or near 



<7 = 0. This criterion is necessary to ensure that all 
free-energy barriers have vanished at this j3, so that the 
different thermodynamic states in the low-temperature 
phase(s) will be visited with a probability corresponding 
to their thermodynamic weight. 

B. Numerical results 

As mentioned in Sec. I, there are indications from 
mean field studies that a SG phase with RSB 

may exist in an intermediate region of the phase dia- 
gram, which roughly coincides with the irreversibility re- 
gion found in dynamical calculations One can see 
a hint of this by doing a Monte Carlo simulation, and 
computing the spectrum of Cy , for a single disorder re- 
alization at a relatively small system size (in order for 
the calculation to be feasible) over the whole phase dia- 
gram. We plot the first two eigenvalues of Cy for a single 
disorder realization in Fig. [|. In the FM and PM phases 
there is clearly a single large eigenvalue Ai, while in a 
region between the two phases the first two eigenvalues 
are of the same order of magnitude. This could be due 
to an intermediate phase with RSB. On the other hand, 
such behavior could also be due to the FM and PM eigen- 
vectors (which represent two distinct forms of ordering) 
trading places in the spectrum as one moves from one 
phase to the other. 

We can follow up these hints by studying larger sys- 
tems, at a single (low) temperature, over a range of fields 
A which spans the intermediate region seen in Fig. |^. 
Hence we have performed a rather thorough Monte-Carlo 
scan of the field parameter A at a fixed, low temperature 
T = 1.5. For a preliminary coarse scan we studied the 
values A/J = 2.0, 3.0, and 4.0. The system sizes for each 
A are L — 5,7,9, and 11, with the number of disorder 
realizations being 2000, 2000, 850, and 650, respectively. 
Previous work |3j has yielded the dynamical phase dia- 
gram shown in the inserts of Figs. ^ and |^. Our coarse 
scan has been chosen to both bracket and sample the ir- 
reversibility region, as shown in the inserts in Figs. ^ 
and ||. In these figures we plot (respectively) the average 
and typical values (where [Ai]t yp = exp([ln Aj] aV e)) for the 
Xi/N, plotted against N on a log scale. There is a clear 
power-law falloff of A2/-/V for A = 3.0 and 4.0. The same 
is not quite true for A = 2.0. That is, for A = 2.0 the 
falloff at small N grows weaker ([A]t yp ) or even ceases 
([A] a „) at larger N. Therefore we performed a second, 
finer scan over the values A = 1.8, 2.2, 2.3, 2.4, and 2.6, 
with the same number of of disorder realizations as given 
previously. 

Results from this scan ([Aj] aV c only) are plotted in Fig. 
||. The hints of possible RSB seen in Fig. [| are consider- 
ably stronger here, particularly at A = 2.4. A = 1.8 is 
clearly in the FM phase. Also, it seems very likely that 
the falloff of A 2 /iV at A = 2.6 will persist at larger N. 
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Hence, we may take A = 2.6 to lie in the PM phase. The 
question is then, have any of the other plots of X2/N in 
Fig. |6| reached asymptotic behavior? 

It is in fact impossible to answer this question from 
examination of Fig. ^ alone. One main problem is the 
rise in X2/N that occurs in many of the plots at small 
N. This obvious finite-size effect is likely due to the fact 
that we are not far from one (or two) phase boundaries. 
Such a rise must be followed by an asymptotic behavior 
which is either flat or falling; and, as is clear from the 
A = 2.6 plot, the latter can occur. For these reasons the 
flat region at A = 2.4 gives no unambiguous conclusion: 
it may be just a local maximum in N, or it may be the 
large- TV behavior. And we are equally unsure of the large- 
N behavior for A = 2.0-2.3. Taken as a whole, the six 
plots in Fig. ^ can be viewed either as an RSB phase 
developing between two other phases, or as the finite- 
size effects of a phase transition between two competing 
order parameters. 

We find that this ambiguity is considerably diminished 
by a study of the equilibrium P(q) distributions for the 
various A values. Figure fj] shows the full P(q) distribu- 
tion for the same A values seen in Fig. H The overall 
impression is clear: two different kinds of order, repre- 
sented by peaks at different values of the overlap q, are 
competing (via fluctuations) in the intermediate values 
of A. At the lower end of the scan, A = 1.8, the FM 
phase gives a peak at q very close to one. There is a 
much smaller peak at q near —1 (not shown in the Fig- 
ure) which comes from fluctuations into the "wrong" spin 
ordering, i.e., that which is not favored by the net ran- 
dom field [which is, again, of 0(N 1 ^ 2 )]. 

At the other extreme (A = 2.6) it is equally clear that 
FM fluctuations (present at small N) disappear at larger 
N, leaving a clear peak in P(q) at the (lower) q value 
appropriate to the PM phase. The dying off of these fer- 
romagnetic fluctuations is well correlated with the decay 
of X 2 /N in Fig. |. 

This interpretation is readily extended to the remain- 
ing A values. For instance, at A = 2.2 in Fig. [?] there 
are clearly PM fluctuations giving a very broad peak 
at q ~ 0.9. This is even more clear if we expand the 
negative-g portion of the P(q) plot for the same A (Fig- 
ure ||). Here we see a small peak near q = — 1 showing 
overlap between "wrong" and "right" FM states, and an 
even smaller peak which we ascribe to overlap between 
"wrong" FM fluctuations and PM fluctuations. Very sim- 
ilar behavior is seen in P(q) for A = 2.3 — both at positive 
and negative q. The main difference is that the PM fluc- 
tuations are stronger than they are at A = 2.2. Hence 
we conclude that both of these A values lie in the FM 
phase for the RFIM, and that the growth of X2/N with 
increasing N for these A values is purely due to finite 
size, plus nearness to a phase boundary. 

Finally, we see no justification for postulating a dis- 
tinct thermodynamic phase based on the behavior at the 



remaining A value (A = 2.4). The behavior both of 
its eigenvalue spectrum, and of its P(q) distribution, is 
clearly intermediate to the behaviors at lower and higher 
A. Furthermore, its P(q) distribution, viewed on its own, 
gives no convincing sign of converging towards a distribu- 
tion appropriate for a spin glass. In fact, the trend with 
N (Fig. ^|) shows that the PM peak is steadily growing 
with TV, at the expense of the FM peak. This is, again, a 
picture of competition between two forms of order (each 
rather strong) near a phase transition, rather than a dis- 
tinct phase. Furthermore, if one plots the total weight at 
negative q [that is, the integral J_ = J dqP(q)] against 
N for the various A values, one finds very similar be- 
havior for A = 2.6 and 2.4: i_ falls very fast, reflecting 
the fact that the "wrong" FM fluctuations vanish rapidly 
with increasing system size. For all of these reasons, we 
tentatively assign A = 2.4 to the PM phase. 

Figs. [| 0, and |[ viewed together, give, we believe, 
a rather clear picture. There are two phases — FM and 
PM— of the bimodal RFIM at T = 1.5 J. The transition 
between these two phases occurs (probably) between A = 
2.3 and 2.4. 

There is a further conclusion which seems obvious from 
the above considerations, and from Fig. ^. The loca- 
tion (in q) of the peak in q, for the FM phase, shows 
no sign of going to zero as the transition is approached 
from below. Hence, either the magnetization is discon- 
tinuous at the phase boundary [fl7|~p2[ or it is continu- 
ous, but with an extraordinarily small exponent j3 Jl5| , ^6[ , 
such that it is 'practically' discontinuous. Obviously, our 
numerical studies cannot distinguish these two possibil- 
ities. However, this aspect of the behavior of P{q) at 
T = 1.5 — that is, the abrupt vanishing of the magnetiza- 
tion at the transition — is completely consistent with the 
behavior seen at higher temperatures (closer to the pure 
critical point) and at zero temperature. This consistency, 
we believe, adds further support to our claim that there 
is no new, third phase in the part of the phase diagram 
studied here. 

IV. CONCLUSION 

The above results and reasoning have presented a 
rather clear picture: there is no spin-glass phase, and no 
replica symmetry breaking, for the bimodal random-field 
Ising model in the low-temperature region that we have 
examined. One of the tests that we have applied was to 
study the scaling with N of the spectrum of eigenvalues of 
the spin-spin correlation function CV, . We have used this 
technique previously |l2] , |l"3| to give unambiguous confir- 
mation of RSB in a case (the Shcrrington-Kirkpatrick 
problem in infinite dimension) where it is strongly be- 
lieved to occur, and also to give strong evidence against 
RSB for a finite-dimensional version of the same problem 
(for which there is still controversy) . Here we have sought 
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evidence for RSB in a different problem, the RFIM — 
which may be mapped to the diluted antiferromagnet in 
an external field, a system which may be studied exper- 
imentally. We have found that the eigenvalue spectrum 
gives no sign at all of RSB over most of the range of ran- 
dom field strength A that we have studied. There is a 
small range of A for which the behavior of this spectrum 
is ambiguous. However examination of the overlap distri- 
bution P(q) removes the ambiguity, giving instead clear 
evidence for a single phase transition between the FM 
and PM phases, with no intervening equilibrium phase. 
Hence we view the behavior of the eigenvalue spectrum 
(Figs. |^ and |^) in this intermediate region as purely the 
finite-size effects of closeness to a phase boundary. In 
particular, the failure of X2/N to decay for this range of 
fields A, even out to N = 1331 spins, stems from the 
persistence of fluctuations between two adjacent phases, 
each with a large order parameter. Hence we speculate 
that the eigenvalues would behave less ambiguously — i.e., 
they would display a reduced tendency to mimic RSB — 
if the transition were continuous in a conventional sense, 
such that the order parameters on either side vanished 
more smoothly. 

It is possible (although, we feel, very unlikely) that 
studies at larger system sizes may reveal these conclu- 
sions to be wrong, confirming instead a spin-glass phase 
confined to a narrow range of A at low T. A more promis- 
ing direction for future work might be to examine the 
RFIM with a Gaussian distribution of fields, since it re- 
mains possible that the two problems have differing phase 
diagrams. Also, one might argue that a SG phase is sup- 
pressed by the strength of the order parameters on each 
side of the transition. By this argument, one should look 
at higher T, closer to the pure critical point. However, 
both nonequilibrium [^|-§} and equilibrium (H| mean-field 
arguments point to the low-temperature region as being 
most likely to support a spin glass phase. This is why 
we looked at low T; and our results rather clearly imply 
only two phases there. 

The authors acknowledge helpful discussions with A. 
Aharony, H. Castillo, and E. Fradkin. GSC thanks the 
Theory Group, Physics Institute, University of Oslo, for 
hospitality. This work was supported by the National 
Science Foundation under Grant DMR-9820816 and the 
Welch Foundation. 



Rev. B 31, 1682 (1985). 
[4] G. S. Grest, C. M. Soukoulis, and K. Levin, Phys. Rev. 

B 33, 7659 (1986); U. Nowak and K. D. Usadel, Phys. 

Rev. B 44, 7426 (1991). 
[5] J. R. L. de Almeida and R. Bruinsma, Phys. Rev. B 35, 

7267 (1987). 

[6] T. Schneider and E. Pytte, Phys. Rev. B 15, 1519 (1977). 
[7] M. Mezard and A. P. Young, Europhys. Lett. 18, 653 
(1992). 

[8] M. Mezard and R. Monasson, Phys. Rev. B 50, 7199 
(1994). 

[9] M. Guagnelli, E. Marinari, and G. Parisi, J. Phys. A 26, 

5675 (1993); D. Luncaster, E. Marinari, and G. Parisi, J. 

Phys. A 28, 4481 (1995). 
[10] E. Brezin a nd C. De Dominicis , Europhys. Lett. 44, 13 

(1998); and |cond-mat/0007457| 
[11] S. Bastea and P.M. Duxbury, Phys. Rev. E 58, 4261 

(1998). 

[12] Jairo Sinova, Geoff Canright, and A. H. MacDonald, 

Phys. Rev. Lett. 85, 2609 (2000). 
[13] Jairo Sinova, Geoff Canright, Horacio E. Castillo, and A. 

H. MacDonald, Phys. Rev. B 63, 104427 (2001) 
[14] See, for ex ample, Fig. 3 of M. Itakura, cond 



[1] For a recent review see T. Nattermann, in Spin Glasses 
and Random Fields, A. P. Young, editor (World Scientific, 
Singapore, 1997). 

[2] H. Yoshizawa and D. P. Belanger, Phys. Rev. B 30, 5220 
(1984). 

[3] C. Ro, G. S. Grest, C. M. Soukoulis, and K. Levin, Phys. 



mat /00 12488 1 which shows the range in the phase dia- 
gram of several previous numerical studies. 
[15] A. Falicov, A.N. Berker, and S.R. McKay, Phys. Rev. B 

51, 8266 (1995). 
[16] A.T. Ogielski, Phys. Rev. Lett. 57, 1251 (1986). 
[17] A.K. Hartmann and U. Nowak, Eur. Phys. J. B 7, 105 

(1999) . 

[18] J.-C. Angles d'Auriac and N. Sourlas, Europhys. Lett. 
39, 473 (1997). 

[19] MR. Swift, A.J. Bray, A. Maritan, M. Cieplak, and JR. 

Banavar, Europhys. Lett. 38, 273 (1997). 
[20] J. Machta, M.E.J. Newman, and L.B. Chayes, Phys. Rev. 

E 62, 8782 (2000). 
[21] H. Rieger, Phys. Rev. B 52, 6659 (1995); H. Rieger and 

A. P. Young, J. Phys. A 26, 5279 (1993). 
[22] A. P. Young and M. Nauenberg, Phys. Rev. Lett. 54, 2429 

(1985). 

[23] S. F. Edwards and P. W. Anderson, J. Phys. F 5, 965 
(1975). 

[24] Analogous behavior occurs in superfluids. Here the quan- 
tity which 'freezes' is a creation or annihilation operator, 
and the analog of our Cy is a one- (bosons) or two- (su- 
perconductors) particle density matrix. The signature of 
the freezing is again an extensive eigenvalue of the den- 
sity matrix. See C. N. Yang, Rev. Mod. Phys. 34, 694 
(1962). 

[25] D.A. Huse and D.S. Fisher, J. Phys. A 20, L997 (1987). 
[26] A. Aharony, Phys. Rev. B 18, 3318 (1978). 
[27] K. Hukushima and K. Nemoto, J. Phys. Soc. Japan 65, 
1604 (1996) ;K. Hukushima, Phys. Rev. E 60, R1008 

(2000) . 

[28] R. N. Bhatt and A. P. Young, Phys. Rev. B 37, 5606 
(1985). 

[29] Helmut G. Katzgraber, Matteo Palassini and A. P. 

Young, |sond-mat/0007113| . 
[30] K. Hukushima, Phys. Rev. E 60, 3606 (1999). 



G 



0.8 



0.6 



0.4 



0.2 



• q „ , eaussian 

. (2) EA 

An gaussian 

(i) 

On bi-modal 

(2) 

A q bi-modal 



100 



1000 



A 

10000 



100000 



t„ 



FIG. 1. q EA and q EA vs. Monte Carlo time to for 
T/J = 2.2, A/J = 0.4, and L = 5. Here we show each 
quantity for the bi-modal (open symbols) and Gaussian (filled 
symbols) random field distribution. 
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FIG. 2. Swapping acceptance ratios among the configura- 
tions at different temperatures at different stages in the itera- 
tion defined by Eqs. (§) and (|). Here L = If, T mi „/J = f .5, 
and A/J = 0.6. 




FIG. 3. The first (upper sheet) and second (lower sheet) 
eigenvalue of the spin-spin correlation matrix for a single dis- 
order realization. The x-axis is temperature (label not indi- 
cated). Here N = f25. The FM phase lies in the front left 
corner (small T and A), where A2 << Ai. The PM phase 
lies on the other side of the "ridge" in A2; in the PM phase 
it is again true that A2 << Ai. The ridge in A2, with the 
corresponding "trough" in Ai, may indicate a distinct ther- 
modynamic phase; or it may simply indicate the effects of the 
FM/PM phase transition. 
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FIG. 4. Scaling of the average of the first six eigenvalues 
\i as a function of N at T/J = 1.5 and A/J = 2.0 (a), 3.0 
(b), and 4.0 (c). The location in the phase diagram obtained 
through dynamical Monte Carlo obtained from reference 3 is 
shown in the insert (here J = 1). 
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FIG. 6. Scaling of the first six eigenvalues (disorder aver- 
aged) at T = 1.5, with A values as shown. The straight line, 
y — l/N, is included as a reference for the eye. 
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FIG. 5. Scaling of the typical ([Ai] typ = exp([ln Ai] avc )) of 
the first six eigenvalues as a function of N at T/J = 1.5 and 
A/J = 2.0 (a), 3.0 (b), and 4.0 (c). The location in the phase 
diagram obtained through dynamical Monte Carlo obtained 
from reference 3 is shown in the insert (here J = 1). Note 
that, in this figure and in Figure the falloff of A2 with N 
becomes weaker at larger N for A = 2.0. 
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FIG. 7. The overlap probability distribution P(q), plot- 
ted for various A (as indicated) and system size. The heavy 
curves correspond to the largest size, L = 11. Curves (light 
lines) for L = 5, 7, and 9 may be located by noting that, al- 
most everywhere in q, the curves monotonically approach the 
L — 11 curve. Note that, for clarity, almost every plot uses a 
different range on each axis. 
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FIG. 8. P(q) for A = 2.2, -1 < q < -0.5. At the largest 
L (L = 11) we see two small (compare the vertical scale with 
Fig. (?]) but distinct peaks, which we ascribe to the overlap of 
"wrong-sign" FM fluctuations with "right-sign" FM ordering 
(leftmost peak), and of "wrong" FM fluctuations with fluctu- 
ations into the paramagnetic phase (peak at q ~ —0.73). 
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